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We numerically solve the transport equations for a quark gas described 
, by the the Nambu-Jona-Lasinio model. The mean field equations of 

■ motion, which consist of the Vlasov equation for the density and the 
00 ' gap equation for the mean field, are discussed, and energy and mo- 

o : 

, mentum conservation are proven. Numerical solutions of the partial 

■ differential equations are obtained by applying finite difference meth- 
CL|| ods. For an expanding fireball of the light quark mass evolves from 
Oh. small values initially to the value of 350 MeV. This leads to a deple- 

' tion of the high energy part of the quark spectrum and an enhancement 

at low momenta. When collisions are included one obtains an equa- 
tion of the Boltzmann type, where the transition amplitudes depend 
on the properties of the medium. These equations are given for fla- 
vor ^[/(S), i.e. including strangeness. They are solved numerically in 
the relaxation time approximation and the time evolution of various 
observables is given. Medium effects in the relaxtion times do not sig- 
nificantly influence the shape of the spectra. The mass of the strange 
quark changes little during the expansion. The strangeness yield and 
the slope temperatures of the final spectra are studied as a function of 
the size of the initial fireball. 
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I. INTRODUCTION 



Deconfinement and chiral symmetry determine the properties of strongly inter- 
acting matter at high temperatures and densities to a large extent. While deconfine- 
ment is yet far from being understood, chiral symmetry restoration and the resulting 
effects can be shown to play the dominant role for the low energy sector of strong 
interactions and can be well studied within phenomenological models. However, 
most of these studies are restricted to the case of thermal equilibrium. Since heavy 
ion experiments are non-equilibrium phenomena, there is a need for incorporating 
chiral symmetry breaking and restoration into a non-equilibrium theory. This is the 
purpose of this paper. 

In order to develop a theoretical model for the evolution of heavy ion collisions 
including a dynamical breaking of chiral symmetry, one has to choose a model 
Lagrangian for the underlying interaction. Here we employ the Nambu-Jona-Lasinio 
(NJL) model in its two and three flavor versions Although nonconfining and 
nonrenormalizable, this model describes well the properties of light hadrons in the 
vacuum. Recent lattice results indicate that also at high temperatures around the 
phase transition the temperature dependence of the condensate and of certain masses 
is correctly reproduced by the model . Its relative simplicity is an additional good 
reason to start from this effective interaction. The NJL Lagrangian is used in the 
Schwinger-Keldysh formalism for non-equilibrium Green functions and equations of 
motion for the particle densities are obtained. Since their derivation is difficult, we 
give only a brief outline of this calculation and refer for more details about this 
subject to Refs. 

In the present work we focus on the numerical solution of the transport equa- 
tions. To this end we firstly investigate the mean field equations of motion, i. e. the 
Vlasov equation, which describes the evolution of the density evolving in a mean 
field, and the gap equation, which describes the selfconsistent coupling of the mean 
field to the density. We show that particle number, energy and momentum are con- 
served. The latter two conservation laws rely on the existence of a potential energy, 
which we derive. The Vlasov equation is a partial differential equation, which we nu- 
merically solve using a finite difference method. For simplicity, we confine ourselves 
to spherically symmetric systems. The numerical solution of the Vlasov equation is 
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presented and the physical effects generated by the mean field are discussed. 

In a second part, binary collisions are included into the equations of motion. 
Since we also want to study the generation of of strange quarks via the processes 
uu — s> ss and dd ss, we use the three flavor version of the NJL model. After 
sketching the derivation of an equation of the Boltzmann type including medium 
effects, we solve this equation using a relaxtion time approach. The relaxation times 
entering this ansatz are computed from cross sections, which depend on temperature 
and are given in Ref. [^]. As for the Vlasov equation, a numerical solution of the 
Boltzmann equation is given. 

We investigate the how the numerical results change when the medium effects are 
switched off in the relaxation times and in the mean field. We investigate also how 
a variation of the size of the initial system influences the final particle spectra and 
multiplicities. In particular, we study the chemical equilibration of strange quarks 
as the initial fireball grows larger. 

To our knowledge, it is the first time that a kinetic equation for a strongly 
interacting system is solved where the input - mean fields and collision terms - 
are derived consistently from an underlying relativistic Lagrangian and where all 
necessary medium effects are included. In particular the underlying symmetry - 
here chiral symmetry - is consistently preserved. Also in the numerical solution we 
have chosen an unconventional way: Solving a partial differential equation instead 
of the test particle method. Of course, also a price has to be paid: The investigated 
system, a quark gas which, at the present stage, does not have the possibility of 
hadronization, is rather unrealistic and any comparison of our results with data 
from heavy ion collisions contains large uncertainties. Therefore the results of the 
paper are important for further theoretical work. 

II. THE VLASOV EQUATION FOR THE TWO FLAVOR NJL MODEL 

In this section, we give a brief derivation of the Vlasov equation and a short 
introduction to the NJL model in its two flavor version and discuss the conservation 
laws. Then we outline the numerical methods used for solving the Vlasov equation 
and present results for one initial condition. 



3 



A. Derivation of the Vlasov Equation 



The starting point of our discussion is the equation of motion for the real time 
Green function G^, which in coordinate space is defined via |^ 

G<{x,y)=i{i,{y)^{x)) , (2.1) 

where iplx) is is the destruction operator for a quark at space-time point x, iply) = 
ip^{y)'y^ and (. . .) denotes the average over the ensemble under consideration. Note 
that is a matrix in Dirac space as well as in color and flavor space. At the mean 
field level, the equation of motion for G^{x, y) can be written as 

{i^^-m,{x))G<{x,y) = Q . (2.2) 

In this equation, mq{x) is a space-time dependent mass, which is given as the sum of 
the bare mass and the Hartree part of the self energy 0,^. We transform Eq. ( p.l|) 
to phase space with the help of the Wigner transformation 

G<{x,p) = J d^e'P'^G^ (^ + ^'^-|) (2-3) 
and introduce the quasiparticle approximation for G'^{x,p) 



G^{x,p) = 2^ (a; p) NNf (^ + "^'?(^)) [^iPo - E{x,p))ng{x,p) (2.4) 

-5ipo + E{x,p)) {2N,Nf - n,-(x, -p)) 



where E{x^p) = ^fP + m'^(x). In Eq. ( ^.4[ ) color and flavor indices are shown 
explicitly and it is assumed that the particle distributions for each individual degree 
of freedom are independent of color, flavor and spin, so that they can be expressed 
by ng{x,p) and nq{x,p), which are the total number of quarks and antiquarks per 
phase space cell, respectively. Inserting Eq. ( |2.4D into Eq. ( |2.1| ) and keeping only 
the leading order in a gradient expansion leads to the Vlasov equation 

dt + v{x,p)d^ - d^E{x,p)dp\ ng{x,p) = , (2.5) 

which is valid for both Ug and rig. Here the velocity is defined by 

P 

v{x,p) = — — - = dpE{x,p) . (2.6) 
E{x,p) 

Note that Eq. (|2.5| ) has been derived without making assumptions about the inter- 
action. It is thus generic for any model H]. 
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B. The Mass of the Constituent Quark in the NJL Model 



In order to be complete, Eq. ( pTSf ) has to be supplemented by an equation for the 
quark mass mq{x), which in turn has to be given by some model for the interaction. 
In the following we will use the SUf{2) version of the NJL model, which is defined 
by the Lagrangian 

,2 / - \ 2" 



(2.7) 



where rrioq is the current quark mass, G a coupling constant and r the Pauli matrices 
in flavor space. For the calculation of the space-time dependent mass mq{x), we 
confine ourselves to a selfconsistent Hartree approximation, which is the lowest order 
in an l/Nf. expansion In this approximation, mg{x) is given by 

mg{x) = mQg — 2iG Tt iC^i^XjX) , (2.8) 

where the trace runs over spin, color and flavor degrees of freedom. After inserting 
Eq. ( p. 41) into Eq. ( p.8| ) and performing the po integration, one obtains the NJL gap 
equation []T|: 

/X , on f \ f "^N^Nf - nq{x,p) - nq{x,p} 
mq{x) = rriQq + 2G'm, x / -— ^ — — . 2.9 

Since the NJL model is not renormalizable, an 0(3) cutoff A has been introduced 
which makes the integral finite. 

C. Conservation Laws 



The coupled equations (|2.5| ) and ( |2.9|) describe the time evolution of the quark 
plasma in phase space. Before we describe the numerical solution of these equations, 
we prove that they conserve the total number of quarks and antiquarks 

AT, = y -J^nqix^p) Mq = j -J^nqix^p) , (2.10) 



the total momentum 



'P= I '^^^^P{nq{x,p) +nq{x,p)) , (2.11) 
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and the total energy S^ot, which will be defined later. 



The time derivative of the total quark number is given, after inserting ( |2.5|) into 

im), by 

^ = / {-vix,p}l + lEix,p)d,) n,{x,p) . (2.12) 

After integrating the second term by parts, the two volume integrals in Eq. ( p.l2|) 
cancel. The remaining surface terms also vanish, if nq{x,p) vanishes both for high 
momenta and large values of \x\. Thus the number of quarks and analogously 
the number of antiquarks is conserved. Note that the conservation of the particle 
numbers can be proven independently of the functional form of the interaction. 

The proof of energy conservation is less trivial and depends crucially on how 
mg{x) is calculated from a dynamical equation. This is obvious since there are 
forms of interactions, e. g. explicitly time dependent external fields, which do not 
conserve energy. In order to prove energy conservation for the Vlasov equation in 
the NJL model, we define the kinetic energy as 

£kin = j '^^^i^E{x,p){nq{x,p) + nq{x,p)) . (2.13) 



After inserting the Vlasov equation (2^), one obtains 
dS]^in f d^x d^p 



dt J {2tt 



dtE{x,p) {ng{x,p) + ng{x,p)) (2.14) 
+ (^-pdx + mq{x)d^mq{x)dp^ {nq{x,p) + nq{x,p)) 



The second term on the right hand side can be transformed to surface integrals and 
does not contribute. The first term can be rewritten to give 

d^kin /■ ,3 ( . [ d^P ng{x,p) +nq{x,p) 

- ' d^xmq[x)dtmq[x) / — — — ^ . (2.15) 



dt J 1^ J - iy ' J (27r)3 E{x,p) 

In order to proceed further, one has to make the crucial assumption that the system 
does not contain quarks with momentum larger than A. Then the momentum 
integral in Eq. ( |2.15| ) can be transformed by using the gap equation (?!?): 

d^p nq{x,p) +nq{x,p) _ f d^p 1 mpg - niq 

(27r)3 E{x,p) ' ^ J\p\<A{2n)^E{x,p) 2Gmq ' ^' ^ 



The right hand side now depends on x only through the mass, so that we can define 
an effective potential Ves via 
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d^p ng{x, p) + ng{x, p) 



dVcsimg{x)) 



(27r)3 E{x,p) ^^{x) 
The effective potential is given explicitly by P,|TO| 



dm„ 



(2.17) 



dfi 



{nig - m^qf 
AG 



2G 



H — - arsinh 



dp 



p 



(-) 



+ ^^ J 

N^NfA 
87r2 



(2.18) 



^A2 + m2(2A2 + ml) 



Note that this result has been already derived in Ref. ||^ from a field theoretical 
point of view. After inserting Eq. ( |2.17| ) into Eq. ( p.l5|) , one obtains 



dt 



d 
"dt 



d xVes{mq{x)) . 



By defining the potential energy E^ot as 



pot 



and the total energy as E^, 



d^xVeE{rnq{x)) 
+ £pou Eq. ( |2.19| ) gives 

O^^tot 



dt 



, 



(2.19) 



(2.20) 



(2.21) 



i.e. the total energy is conserved. It can also be read off from Eq. ( 2.19| ), that the 
quantity 

d^p 



t{x) 



;E{x,p) {ng{x,p) + ng{x,p)) 



+ Vesimqix)) 



(2.22) 



has the meaning of an energy density. We note again, that Eq. ( 2.21 ) can only 
be obtained by assuming that no particles with momentum larger than A appear 
in the system. This somewhat artificial restriction arises as a consequence of the 
introduction of the 0(3) cutoff into the NJL model. The form of Ves as a function of 
the constituent quark mass rrig, as evaluated from Eq. ( p.l8| ), is shown in Fig. |l]. Vcs 
has a minimum at the constituent quark mass in the vacuum, since, according to 
Eq. ( |2.17|) , VcfT arises as an integral over the vacuum gap equation. In the vacuum, 
the gap equation is thus equivalent to the condition that the energy density, as 
defined by Eq. ( |2.22| ), is minimized. In Fig. |I], we have also added a constant to 
Eq. (glSD , so that Ves = at the vacuum quark mass and thus vanishes in the 
vacuum. 
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As for the case of energy conservation, momentum conservation can also be 
proven only if the interaction is specified. The time derivative of the total momentum 
is obtained by inserting (|2.5| ) into After performing an integration by parts 

and dropping the surface terms, one obtains 

'^ = - j'^0^{n,{x,p)+n,{x,p))d^E{x,p) . (2.23) 

The structure of the momentum integral is now the same as in Eq. (|2.15| ). One thus 
has 



Itt 



J (fxd^Vesim,{x)) = , (2.24) 



i. e. the total momentum is conserved. As Eq. ( p.21| ), Eq. ( p.24|) is only valid strictly 



if no particles with momentum larger than A are present. 



D. Numerical Methods 



The Vlasov equation (|2.5| ) is a first order partial differential equation, which 
in general depends on the seven variables t, x and p. The usual solution method 



for this equation is the test particle method |[TT|,[T2|, which is based on a solution 
of the characteristic equations. Here we do not follow this way, but rather solve 
the equation by a finite difference method in order to study the advantages and 
disadvantages of this approach. In order to reduce the numerical effort, we assume 
that the system is spherically symmetric. This, in turn, implies that ng{x,p) can 
only depend on time and the three variables 

xp 

r = \x\, P=\p\ , V = T^TT^ ■ (2-25) 

fI In 

The last of these variables, rj, is the cosine of the angle enclosed between x and p. 
With the variables of Eq. (|2.25|) , Eq. ( |2.5| ) transforms to 



dtrig = ^ 



( » \ > A I ^qdrmq p\ 
T] [rnqdrmq dpUg -pdrUq) + [I - Tj j I -I OrjUq 



(2.26) 



where the arguments of ?7,g, and E have been dropped for simplicity. In order 
to solve Eq. ( ^.261 ), we use a finite difference scheme. In detail this proceeds as 
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follows: Firstly, the variables r, p, rj are discretized equidistantly according to the 
prescription 



ri = ir^e.jNr , i = 1 . . . Nr (2.27a) 

Pj = J Pm.jNp , j = l...Np (2.27b) 

r^fc = [2(A: - l)/(iV, - 1)] - 1 , k = 1 . . . . (2.27c) 

This parametrization avoids the coordinate singularities of Eq. ( p.26|) occurring at 
r = and p = 0. The time derivative of Uq at time tn is replaced by 

dtnq{t,r,p,r]) ^ ^[nq{tn+i,ri,pj,r]k) - ng{tn,ri,pj,r]k)] . (2.28) 

In order to achieve numerical stability, we use a fully implicit procedure for the 
representation of drUg at time tn, 

{■^[nq{tn+i,ri+i,pj,rik) -nq{tn+i,ri,pj,r]k)] , i = l 
j^^[nq{tn+i,ri+i,Pj,r]k) -nq{t„,+i,ri_i,pj,r]k)] , 1 < i < Nr 
■^[nq{tn+i,ri,pj,r]k) -nq{tn+i,ri_i,pj,r]k)] , i = Nr 

(2.29) 

and analogous prescriptions for dpUq and dr^Uq. Note that this implies open boundary 
conditions, i. e. we assume that Uq can be continued smoothly to points outside 
the grid. The time update for Uq is done using an operator splitting scheme |]13 



i. e. we sequentially perform the three update steps which would arise if only the 
terms proportional to drUq, dpUq and i9^nq were present on the right hand side of 
Eq. ( p.26|) . Technically, each update consists of the solution of tridiagonal linear 
systems, which can be implemented efficiently using the Thomas algorithm []T3|. 
With these ingredients at hand, the calculation proceeds as follows: 

i. Initialize Uq and riq according to some given initial conditions. 

ii. Compute niq at each space point. 

iii. Compute physically relevant quantities like particle density, energy density 
etc. and store them for later evaluation. 

iv. Calculate Uq and riq on the next time slice. 
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V. Proceed with step p. 

The calculation stops at that time, when the system has become a system of nonin- 
teracting particles. A sufficient criterion for this is that the quark mass comes close 
to its vacuum value everywhere. 



E. Numerical Results 



In this section we discuss our numerical results for the SUf{2) Vlasov equation. 
The model parameters chosen are rriQg = 5.0 MeV, A = 653 MeV and GA^ = 2.10 
T^. The initial conditions correspond to a spherically symmetric ffieball around 



the origin with baryon number zero: 

nq{x,p) = nq{x,p) = exp -—^ /, 



(2.30) 



with a cutoff function fc{x) = (1 — tanha:;)/2. The ffist factor of Eq. ( 2.30| ) describes 



the spatial shape of the system. The second factor is a Fermi distribution, which 
describes the distribution in momentum space. The mass is chosen to be the 
thermal mass corresponding to Tq, i. e. we have complete thermal equilibrium in 
the centre of the fireball. The third factor serves to cut off high momenta. This 
is necessary, since the boundary condition of Eq. ( p.29|) assumes that the particle 
distribution can be continued smoothly beyond the grid boundary. In order to avoid 
the cutoff artifacts which were discussed in Section |11 (J| , and to obtain exact energy 
conservation, it is thus necessary to cut off the distribution smoothly. The function 
fc is designed to interpolate smoothly between the values fc[{p + Pc — -A)/6p] ~ 1 for 
p <^ A — pc and fc[{p + Pc — A.) / 6p] = for p ^ A — pc, the drop between these two 
values taking place within an interval 6p around p = A—pc- The standard parameters 
chosen for Eq. (|2.30|) are ro = 3 fm for the initial radius and Tq = 240 MeV for 



the initial temperature. The parameters entering the momentum cutoff factor are 
Pc = 100 MeV, Sp = 20 MeV. 

The maximal values of r and p on the grid are rmax = 10 fm, Pmax = A. The size 
of the grid is Nr = Np = 100, = 50. 
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The fireball described at t = by Eq. ( |2.30 ) expands for t > due to thermal 
motion. After a certain time the density of particles is sufficiently low that the 
system behaves as a system of non interacting particles. This is demonstrated in 
Fig. ^ where the quark density at various times is shown as a function of r. At 
t = 0, the density is described by the factor exp (— r^/2rg) in Eq. ( ^.301 ). At later 
times, the system expands and the density drops. Although Fig. |^ shows the density 
up to t = 8 fm only, the calculation has been extended to t = 17 fm/c, when the 
density practically vanishes. 

The quark mass mq{r), shown in Fig. |^, behaves in a similar way. At t = the 
quark mass forms a "potential well" due to the high density in the centre. As the 
density drops, the quark mass approaches its vacuum value of m^ vac = 312 MeV. 
For t > 17 fm/c, mq{r) is practically constant. This leads to a vanishing mean field 
term in Eq. ( |2.5| ) and the system decouples. 

The final "observed" momentum spectra are calculated using the observation 
that there exist closed surfaces in coordinate space with the properties that (i) at 
t = all particles are contained within this surface, (ii) at time t = oo all particles 
are outside the surface and (iii) outside the surface the system behaves effectively 
as a system of non-interacting particles. The particle spectra then can be obtained 
by integrating the fiux through such a surface. Choosing the surface to be a sphere 
with radius ri, one obtains 



d^^g 



CO p 

dt / dnnqv± , (2.31) 



d^p (27r) 

where v± is the velocity component perpendicular to the surface. Since in practice 
the time integration in Eq. ( p.31|) can only be extended to finite times, we add 
also the (small) contribution of those particles, which still remain inside the surface 
r = ri when the calculation stops. In Fig. ^, the momentum spectra of quarks are 
shown for the initial and final states. At t = 0, the spectra are governed by the Fermi 
part of Eq. (^^301) up to p ~ 550 MeV, where the distribution is cut off by the factor 
fc [{p + Pc ~ ^)/^p] in Eq. ( p.30| ). During the evolution, the number of particles with 
low momenta is enhanced, while the high momentum region is depleted, since the 
generation of a dynamical mass during the evolution slows particles down and thus 
enhances the low momentum region of the spectrum. 

The time development of the total energy is interesting in two respects: firstly 



11 



as an accuracy check to the program and secondly for a study of the importance 
of the potential energy, which is usually neglected. The time behaviour of £^kin and 
£tot is given in Fig. ^. As was detailled in Section |II Q , the potential energy is 
positive and vanishes in the vacuum. This behaviour can be observed in Fig. ^: 
At small times, when the system consists of low mass quarks, the potential energy 
is maximal and contributes about 5% to the total energy. During the expansion, 
the quark mass approaches its vacuum value and the the potential energy goes to 
zero. Concomitantly, the kinetic energy rises. Calculating kinetic and potential 
energies numerically, we find that energy is conserved to an accuracy of 1%. We 
also obtain particle number conservation with an accuracy of 6%. Both results point 
to a numerical accuracy of the order of 5-10%. The total momentum, on the other 
hand, vanishes identically due to the spherical geometry. 



III. THE BOLTZMANN EQUATION FOR THE THREE FLAVOR NJL 

MODEL 

The results of the previous section are generalized in two respects: by including 
strange quarks and by going beyond the mean field level and adding the collision 
term. The exact expression involving collision integrals is replaced by the relaxation 
time approximation and numerical results are obtained. 



A. Derivation 



After studying a collisionless system, we extend our approach in order to include 
collisions. These can be incorporated by generalizing Eq. (|2.2|) to 0] 



{i^^-m{x))G<{x,z) = j dS {T.\x,y)G<{y,z)-i:<{x,y)G\y,z)) , (3.1) 

where the causal self energy S^, the anticausal Green function and the self energy 
appear on the right hand side. Since we also consider the production of strange 
particles, we have to extend the underlying dynamical model. The Lagrangian for 
the SUfiK) NJL model reads 



f=u,d,: 



a=0 



(3.2) 



K [det ^(1 + 75)^ + det ^{l - 75)^ 
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Equation ( p.2|) contains two coupling constants G and K. The matrices A° are 



the Gell-Mann matrices in flavor space with A° = ■y/2/3. The Green function G< 
in the quasiparticle approximation is still diagonal in flavor space, but no longer 
independent of flavor: 

G^ix,p) = ^^^^^Jj^(^ + mf{x)) [6{po-Ef{x,p))nf{x,p) (3.3) 



-6{po + Ef{x,p)) {2N^ - nj{x, -p) ^ 

where nf{x,p) is the density of quarks of flavor / for / = u, d, s and Ef{x,p) = 
fp' + m?j{x). The mean field derived from the Lagrangian (|3.2|) becomes ||l|JT5|| 

mf[x) = niof + AGmf[x) ——- — — (3.4) 

W|p1<a (27r)3 Ef{x,p) 



d?p 2Nc — nf>{x,p) — np{x,p) 
p1<a (27r)3 Ef{x,p} 

Sp 2Nc — n jn (x, p) — njii {x, p) 



+2Krafi{x) 

'''"^"^'^^Vl<A(2vr)3 Ef„{x,p) 

where /, /' and /" are pairwise distinct flavors. 

After performing a Wigner transformation and a gradient expansion of Eq. ([3] 
inserting Eq. ( |3.3D and taking the real part of the resulting equation, one arrives at 
a Boltzmann like equation [Q 

(dt + Vfd,-d,Efd,)nf{x,p)= {lTh^f'f[+I}f-.f'fi) ' (3-5) 

fiJ'Ji 

where lYfl^ff^ and IJjl^fif are the collision integrals due to quark-quark scattering 
//i — > f f[ and quark-antiquark scattering //i f f[, respectively. One has 

I}f^f'fi = ^ / dQ{2nr6\p + p,-p'- p[)l\Wff,-.f'f^ (3.6) 
nf {x, if)nf^ (x, i^i)(f)f{x, p)(f)f, (x, pi)-nf{x, p)nf^ (x, pi)(j)f (x, p)(j)f^ (x, p'^) 

with the statistical factor 1/2 appearing in front of the squared transition amplitude. 
In Eq. ( |3.6| ) we have used the invariant integration volume 

= ci^Pi d^p' d^p'i . 
^ (2vr)32i?^, {2nr2Ef, {2ny2Ej. ^ ' ^ 

and the blocking factors 
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M^^P) = ^-'%^ ■ (3.8) 

In contrast to Eq. (|3.6| ), the collision integral due to quark-ant iquark scattering does 
not contain the statistical factor. One has 

Iffl^rn = ^ / dQi^-mP + PI-P'- Pi)\Wff-.^f'fi (3.9) 



Note that, in terms of the differential cross section, one has [0 

^ -dQ{2n)^5\p + p,-p'~p[)K\M^=^^dnv,,,^ , (3.10) 



2Ef ' ' ' " ' ' (27r)3 

where k is a statistical factor which equals 1/2 for the elastic scattering of identical 



particles and 1 for all other processes. Equations (|375|), ( p.6|) and ( p.9|) form the 
basis of our following investigations. 

B. Relaxation Time Approximation 



The collision integrals in Eq. (|3.5| ) contain the differential cross sections do /dVt 
for elastic scattering processes like uu uu and inelastic ones like uu — ss. The 
scattering amplitude M. in Eqs. (p.6|), ( |3.9| ) depends on the selfconsistent masses 
mf{x) and on the actual densities nf{x,p) via the Pauli factors 1^,^. Thus has 



to be recalculated at every space-time point of the expansion. This is a formidable 
numerical task and goes beyond the scope of this paper. Therefore we use the 
relaxation time approximation. Our ansatz for the transport equation is 



{d, + v,d^ - d^E.d,) n,{x,p) = ^-^l^^^) y^P) _ IhflR + G,{x,p) (3.11a) 

[dt + Vsd^ - d^Esdp)ns{x,p) = — — - + Gs{x,p) . (3.11b) 

As in Section m, we set nu{x,p) = nd{x,p) = nq{x,p)/2 and also nf{x,p) = nj{x,p). 
The first term on the right hand side of Eq. ( |3.11a|) describes gain and loss terms 
due to the elastic scattering of light quarks. In this term, hq{x,p) denotes an ef- 
fective equilibrium distribution function to be detailled below. The second term 
describes the loss of light quarks due to the processes uu ss and dd ss. The 
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last term models the gain of light quarks from strange quarks due to the reverse 
processes. Its concrete form will also be given below. The terms on the right hand 
side of Eq. (|3.11b| ) have similar origins, with light and strange quarks being suit- 
ably exchanged. The relaxation times r^^, r^s, Tgg and Tgg are calculated from the 
equilibrium cross sections [Q. They will be given in Section [ill (J| . 

The effective densities nq{x,p) and ns{x,p) are chosen to have the form 

"^^^'^'^-exp(pX(^)/T(a:)) + l ^^"^^^^ 

2N 

where u'^{x) is the collective velocity of the fluid component with flavor / at space- 
time point X, T{x) an effective flavor independent temperature and ({x) a factor, 
which measures the amount of chemical equilibrium for the strange quarks. The 
velocity can be calculated from the numerical solution nf{x,p) at every time 

slice, 

"'^^^ = / (^^^^^^'^"^ / / (^^^(^'^"^ ' ^^-^^^ 

so that T{x) and ({x) remain to be determined. 

The requirement of energy conservation is used to fix the effective temperature 
T{x) and the density factor ({x). After introducing a potential energy which, as 
for the SUf{2) case, compensates the change of the kinetic energy due to the mean 
field, one obtains 

^ dStot f d^x d'^p Eq{x,p} , ^ / ^\ T A\ 

{ng{x,p) - nq{x,p}) (3.14) 



2 dt J (27r)3 Tqg{x,p) 

f d^xd^pEs{x,p} , ^ I 
+ / T^TTs 1 — ^{ns{x,p) -n,{x,p)) 

d^xd^p ( Eg{x,p)ng{x,p) 



(27r)3 \^ ' ' rqs{x,p) 

d^xd^p ( „ , ^\r-i f -\ Es{x,p)ns{x,p) 



f a^x d p / ^ , ^ / ^ 



(27r)3 ^ .V ^^^^^^^ J 

for the change of the total energy due to the coUisional part of Eqs. ( 3.11|) . The 
factor 1/2 at the left hand side accounts for the antiparticle degree of freedom. Note 
that for the Lagrangian ( |3.2| ) an effective potential, as it was given in Eq. ( |2.18| ) for 
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the SU{2) case, exists in a mathematical sense, however, due to the six fermion 
couphngs in the 't Hooft determinant, it is not possible to compute it analytically. 
For the exact Boltzmann equation ( p.5|) , each of the integrals on the right hand side 
of Eq. ( p.l4| ) vanishes individually, even if the spatial integrations are dropped. We 



adopt this stronger constraint also for Eq. ( p.l4[ ), and have the two conditions 

{ng{x,p) - ng{x,p)) = (3.15a) 



(271)^ Tyq{x,p) 



d^p E,{x,p) ^Qie;M 

[ns[x,p) -ns[x,p)) =[) (3.15b) 



{2'kY Tss{x,p) 

to determine T{x) and C(x) at each space-time point. In order to fulfill the remaining 



constraints 



(E^{x,ms{x,p)-^^^^^^P^\ =0 (3.16a) 



(27r)3 \ ' ' Tgs{x,p) 

f d^p ( Es{x,p)ns{x,p)\ /Qifiu\ 

J (M^ [E.i^^P)G,{x,p) - j = , (3.16b) 

we make the ansatz that the gain terms Gq and Gs are proportional to the densities 
in local thermal equilibrium: 

G,(x,p)=7.(x)^^4^ (3.17a) 

G,{x,p) = 7.(x)^^^^lMM . (3.17b) 
rsq[x,p) 

The factors 7g(a;) and 7s(x) are determined by solving Eqs. (|3.16|) : 

[ d^p Es{x,p}ns{x,p} I f d^p Eq{x,p}nq{x,p} 

^'^^^-J{2nr r,,{x,p) / J {2nr rq,{x,p) ^'-''^^ 
d^p Eq{x,p)nq{x,p) I f d^p Es{x,p)ns{x,p)/C{x) 



' J (27r)3 Tq,{x,p} I J (27r)3 r,,(x,p) ^ ^ 

With this form of the gain terms, the collision term vanishes in local thermal equi- 
librium. 



C. Relaxation Times 

The basic assumption of the relaxation time approximation is that the system 
is close to equilibrium and that the right hand side of the Boltzmann equation is 
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only linear in the difference nj — nj. Therefore the relaxation times have to be 
evaluated with the equilibrium densities and masses and depend only on the local 
temperature. 

The collision term on the left hand side of Eq. ( p.5| ) involves a sum over all pos- 
sible scattering processes. The task of computing the relaxation times for Eq. ( ^.111 ) 
is thus twofold: firstly one has to define the relaxation time for one single process 
and afterwards add up the contributions of several processes in order to obtain 



Eqs. (KU 



In order to define the the relaxation times, we go back to Eq. ( |3.5| ) and write it 
in the form 

{d, + vjd, - d^Efd,) nf{x,p) = Y: (4,gain - 4,ioss) , (3-19) 

p 

where the sum runs over all possible processes and /pgain ^p,\oss ^^^^ respective 
gain and loss terms due to process P. We focus on one generic process P : ffi ^ 
f'f[. The loss term due to this process is given by 

W) 

By comparing this to the ansatz one immediately obtains the relaxation time 

for this process 

1 f d^pi dap 



Mx,p) J (2vr)3^^"-^^"^^'^^"^^^ • (3.21) 
In order to evaluate this expression, we go to the rest frame of the plasma and 
approximate the differential cross section by its equilibrium value at the local tem- 
perature T(x), the density of flavor /i by 

exp [Ef^{T{x),pi)/T{x)\ + 1 

where Ef-^(T{x),pi) is computed using the equilibrium mass corresponding to the 
effective temperature T{x). The blocking factors are replaced by 

0/.(T(x),s) = 1 f i , (3.23) 

exp[E'f{T{x),s)/T{x)]+l 

where i?j™(T(x), s) is the energy for the participant /' in the centre of mass frame, 
which can be expressed as a function of the masses and the Mandelstam variable s. 
These replacements lead to 
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^ Vrciap^e{i{T{x),s)nf^{T{x),pi) . (3.24) 



TpiTix),pj) J (27r)3 
The effective cross section in Eq. ( |3.24] ) is defined by 

ap,efT(r(x),s) = 0/Kr(x),s)0y,(T(x),s) J dQ^{T{x),s,t) . (3.25) 



By introducing the Mandelstam variable s as integration variable, Eq. ( |3.24|) can be 
reduced to 

1 1 1 ,oo 



ds,J[s - (m/ + m/J2][s _ (^^ _ ^^J2] 



xap^ca(T{x),s)ujff^(T{x),s,Ef) , (3.26) 

where 

^//i(^(2;), E/) = J dEf^nf^{T{x),pi)Q {2pfPf,f - (s - mj - m\ - 2EfEf^)^ 

(3.27) 



is a weight function, which can be computed analytically for rij^ given by Eq. ( p.22| ) 



The temperature dependent cross sections entering Eq. ( p.26| ) are taken from Ref. |Q. 



Note that Tp in the rest system of the plasma is a function of T(x) and pf only. It 
can thus be easily obtained via table lookup. In the lab system, we compute Tp{x,p) 
by applying a Lorentz boost to Eq. ( ^.2(j| ). 

In order to relate the relaxation time for a specific process to those needed in 
Eq. ( p.ll|) , we consider again the loss term of the exact equation (|3.5| ). To be specific, 
we set f = u. Then the exact loss term can be decomposed into the contributions 

P.loss tin^nu, loss ' ud^ud,loss ' us^us,loss yo.^Oj 

P 



= r 

uu- 


1 TU 

-^nUjloss ' ud- 


-> lid, loss 


+ r 

' us- 


-+iiS,loss 


+ r- 

' uu- 


_|_ r"_ 

tin, loss ' ud- 


-+ lid, loss 


+ r- 

' us- 


->MS,loSS 


+ 

' uu- 


1 r« 

■^ddyloss ' uu- 


->ss,loss 







Only the last two of these contributions correspond to inelastic processes. However, 
we have considered u and d to be degenerate in Eq. (|3.11|) , so that we have to include 
also the contribution of uu dd to the elastic relaxation time Tqq of Eq. ( p.ll|) . Thus 
one has 

11 1 C 1 1 C 1 . s 

— = + + + + — + + : (3.29a) 

^^qq ^uu—*uu ^ud—*ud Tus-^us Tuu^uu ^ud^ud Tus^us ^uu^dd 

— = , (3.29b) 
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where a factor ^ has been multiphed to the contributions of the processes us — > us 
and us -h> us in order to relate Eq. (|3.22| ) to Eq. ( |3.12b| ). By considering the loss 
term for the strange quark, one obtains in an analogous fashion 

— = — ^ + + — ^ + (3.30a) 

1 2C 



(3.30b) 

where the factor 2 in the contributions of su — su, su su and ss — >■ uu accounts 
for the fact that in all these processes u can be replaced by d. 

In Fig. ^, we show the resulting inverse relaxation times for the elastic scat- 
tering of light quarks (solid line), the elastic scattering of strange quarks (dashed 
line) and for the production of strange quarks (dotted line) as a function of tem- 
perature and for a momentum p = 200 MeV. While t~^^ and are approximately 
equal, is around two orders of magnitude lower than for this value of p. 
The main reason for this is the threshold for the process qq ss, which makes 
this process improbable at low momenta. The ratio between and r^"^ decreases 
for increasing momenta. At low temperatures, goes to zero rapidly, mainly be- 
cause of the density factor in Eq. ( |3.24| ). The solid, dashed and dotted curves of 
Fig. H are calculated using temperature dependent quark masses and the tempera- 
ture dependent cross sections given in Ref. 0. These are to be contrasted with the 
dot-dashed line, which gives r^"^ from a calculation, where the masses have been 
set to the current quark masses and the cross sections have been calculated in the 
Born approximation [Tl|, which does not include any medium dependence. At high 
temperatures, this procedure results in an inverse relaxation time, which is approx- 
imately a factor 3.5 smaller, whereas the falloff at low temperatures is weaker due 
to the lower masses. 

In Fig. 1^ Tg^^ is given as a function of temperature for p = (solid line) and 
p = 200 MeV (dashed line). One notes that the solid line displays a pronounced 
maximum at T ^ 230 MeV. It originates from a divergence of the scattering length, 
which occurs for the quark-ant iquark scattering processes at the Mott temperature, 
where the pion becomes unbound P,p!7|. This phenomenon is related to critical 
scattering. The Mott temperature is ca. 210 MeV for our parameter set. As to be 
expected, the Born calculation for p = 0, which is given by the dotted line of Fig. |^, 
displays no structure at this temperature. 
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D. Numerical Results 



1. Densities, Quark Masses, Spectra and Strangeness Yield 
The transport equations are solved using an obvious generalization of the 



numerical methods outlined in Section [II D| . In this subsection we show results of a 
calculation with the initial condition 

n,(a;,p) = exp ( I ^ ^ ^ (3.31a) 



0/ e^^[^Jp^ + rni/Toj + 1 
n,(x,p) = . (3.31b) 



Comparing the form of Eq. ( |3.31a| ) with Eq. ( ^.3(J| ) there is only one difference 



namely the absence of a cutoff in momentum space: While Eq. (|2.3CI|) contains a 
smooth cutoff, we now continue the Fermi shape of the distribution up to the grid 
boundary Pmax, which in our standard parameter set is set to Pmax = A. This 
modification is necessary if one wants to determine effective temperatures from an 
exponential fit to the particle spectra. The disadvantage of this modification is that 
the total particle number and the total energy are no longer exactly conserved. We 
will discuss the magnitude of these effects later. 

The standard values of the parameters are the same as used in Section p.lE: 



ro = 3 fm and Tq = 240 MeV; the initialization mass mi is taken to be a solution 
of Eq. ( p.4|) at r = 0. The NJL model parameters for all SUf{3) calculations are 
taken from Ref. pl and are mog = 5.5 MeV, = 140.7 MeV, GA^ = 1.835, 
KA^ = 12.36 and A = 602.3 MeV. 

The time dependence of the light quark density is shown in Fig. ^. This figure 
should be compared with the same quantity calculated for the Vlasov equation, 
shown in Fig. ^ The time steps coincide in both figures. The behaviour of the 
light quark densities is rather similar. The time variation for the density of strange 
quarks is shown in Fig. |^. When comparing it to Fig. ^ for the light quarks, one 
should first note the different scale on the ordinate (about a factor of 30 down), and 
secondly the initial condition Usit = 0,x, p) = 0. The density rises due to strange 
quark production until t = 4 fm/c, after which time the rarefaction due to flow 
overcomes the production and the density drops. 
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The light and strange quark masses are shown in Fig. 10 as a function of time. For 
the hght quark mass one obtains essentially the same behaviour as observed in the 
SU (2) case shown in Fig. ^. At t = the light quark mass amounts to rriq = 58 MeV 
in the centre and rises to the vacuum value nig = 368 MeV in the outer regions. 
As time progresses, the masses increase everywhere to this value. For the strange 
quarks, on the other hand, one does not observe strong variations of the constituent 
mass. This different behaviour in the space-time behavior of mq{x, t) and ms{x, t) is 
significant. While the difference in the vacuum values uis — mq = 182 MeV compares 
well with the mass difference between lambda and proton, m/^ — nip = 177 MeV, 
the in medium corrections in the plasma are large for the light quarks and small 
for the strange quarks. In the NJL model the medium effects on the masses can be 
traced to the Pauli-blocking, essentially. The Pauli-blo eking for the strange quarks 
is small because of the low density of this kind of quarks. At t = 0, the strange quark 
density vanishes and only the term proportional to K in Eq. ( p. 41 ) gives corrections 
to the vacuum mass. This leads to a value of nig = 465 MeV in the centre and 
rris = 550 MeV in the outer regions. At later times, the term proportional to G 
also gives medium corrections. However, these are not sufficient to lower the strange 
quark mass substantially in the central region. With increasing time, also the strange 
quark mass evolves to its vacuum expectation value everywhere, as expected. The 
large mass difference of the in medium masses suppresses the production of strange 
quarks in our calculation - and possibly in nature. 

In Fig. |Tl|, we show the "observed" (for t oo) distributions of particles in mo- 
mentum space, dM/(f J)^ clS db function of the particle energy. An exponential fit is 
applied to the high energy part of the spectra (dashed lines). The effective slope tem- 
peratures of these fits for the light and strange quarks are different: Tg = 185 MeV 
for light quarks and Tg = 164 MeV for the strange ones. The calculated spectra 
show an enhancement at low energies compared to the exponential fit. This en- 



hancement arises from the mean field, as has already been explained in Section [II E . 
We comment on the difference Tg — Tg in Section 

We find for the strangeness yield, defined as the ratio of the number of strange 
quarks to that of light quarks at t = oo, a value of J\fs,t=oo/J^g,t=oo = 2.2%. This 
value is rather low and is mainly determined by incomplete chemical equilibration, 
as we will show in in Section 
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2. Dependencies on the Grid Settings and Numerical Accuracy 



In order to see the reliability of the results and the stability with respect to the 
numerical scheme, we have investigated the dependence of the particle multiplicities 
and the spectra on the grid extension. The results of the various runs are compiled 
in Table |I[ The calculated multiplicities of the light and strange quarks with the 
standard parameters described above are given in the row labeled (a). One has 
Afq^t=o = 343.2 in the initial state and Afq^t=oo = 365.6 and A/'s,t=oo = 8.1 in the final 
state. This corresponds to a relative yield of J^s,t=oo/-^q,t=oo = 2.2%. We observe 
an increase of the total multiplicity from {Afg + ■Afs)t=Q = 343.2 to {Afg + Afs)^^^ = 
373.7. However, A/g + A/'s should be a constant of motion for the exact Boltzmann 
equation (|3.5| ). 

The possible reasons for the increase of the increase are the following: (i) numer- 
ical inaccuracies, (ii) the relaxation time ansatz, which, however, should lead to a 
decrease of the particle number due to the dropping temperature and (iii) the exten- 
sion of the grid in momentum space. The last effect can be understood as follows: 
In the calculation related to row (a) of Table |, the grid has maximal momentum 
Pmax = A = 602 MeV and the initial "thermal" distribution of Eqs. ( p.31|) is cut 
off at this value. Calculating also the particles from Eq. (|3.31| ) with p > Pma.x, we 
see that our grid includes only 40% of all particles explicitly. Nevertheless, due to 
the open boundary conditions which are chosen in Eq. ( |2.29 ), i.e. smoothness at 
P = Pmax, the numerical procedure assumes that the particle density can be extrap- 
olated analytically beyond the grid boundary and thus "recognizes" the remaining 
particles implicitly, at least as far as the derivative terms of the transport equation 
are concerned. Since the mean field tends to slow down particles, it is possible that 
particles with p > Pmax, which are not counted at t = 0, appear inside the grid at 
later times, and thus the multiplicity increases. 

In order to eliminate the influence of the finite grid extension in momentum space, 
two other settings have been used: firstly, the grid is extended to include momenta up 
to Pmax = 3A. In this calculation, all momentum integrals have been extended top = 
Pmax5 except those occurring in the gap equation (p^). The resulting multiplicities 
are given in row (b) of Table |. Compared to the 8.9% increase observed in the 
previous calculation, one now obtains a increase of 5.8% of the particle number. 
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which is comparable with the Vlasov result. The fraction of strange particles in this 
calculation is 3.9%, compared to 2.2% for calculation (a). 

In a second investigation a smooth cutoff for the high momentum region is in- 
troduced in the initial conditions, as is done for the Vlasov equation in Eq. ( p.3(]| ). 
In this calculation, the effective distribution functions of Eqs. ( ^.12| ) are consistently 
cut off by the same factor. This procedure is also consistent with the requirement 
of energy conservation, as has been explained in Section |II C[ and with the nature 
of the NJL model as a low energy theory. The multiplicities for this calculation are 
given in row (c) of Table |. We observe a 1.3% decrease of the particle number, 
which means that the particle number is in good approximation conserved. The 
strangeness content in the final state amounts to 2.3%. 

The light quark spectra of calculations (a), (b) and (c) are shown in Fig. |T2|. The 
data for the calculation (a) (solid line) agree with those for calculation (b) (dashed 
line). The data for the calculation (c) are given by the dotted line. Although 
calculation (c) agrees reasonably with (a) and (b) in the low momentum region (as 
it should), the momentum cutoff factor distorts the high momentum part completely. 
It is thus impossible to fit a reasonable slope temperature to this curve. 

While the agreement of the spectra is good for light quarks, this is not the case 
for strange quarks, as can be seen from Fig. |T3[ Calculation (b), i.e. Pmax = 3A, 
gives a twice higher value for dAfg/d^p as compared to calculation (a) and (c), i.e. 
Pmax = A, since more light quarks with high energy are available for the production 
of strange quarks. 

This may be the appropriate place to report the numerical effort and accuracy. 
The use of a partial differential equation in four variables for the expansion of a 
spherical system in phase space is limited by the use of computer memory and CPU 
time. For a 50^ grid, our program uses ca. 3 MByte memory and 20 seconds of CPU 
time per time slice on a Sparc 20 workstation. The grid size as well as the number 
of time slices has to be increased when e. g. increasing the radius tq of the fireball, 
so that the above numbers have to be scaled accordingly. Most of the CPU time is 
consumed by the determination of the effective temperature. All indicators point to 
an overall numerical accuracy of a few percent of our results. 
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3. Dependencies on Medium Effects 



The inclusion of medium effects in the colhsion terms is an important and new 
aspect of our calculation. Indeed, the medium effects are large, both for the relax- 
ation times shown in Fig. 0, as well as for the quark masses shown in Fig. |10|. It is 
thus important to investigate their influence on the final spectra. Our results are 
compiled in Table y, where row (a) relates to the standard calculation. The medium 
effects are studied in two ways: by keeping or switching off the mean field and by 
replacing the cross sections in the relaxation times by their Born approximation 
[I^ values. All other parameters correspond to the standard set. Multiplicities and 
slope temperatures for these runs are given in rows (d) and (e) of Table y, respec- 
tively. While the results of calculation (d) are similar to those of calculation (a), 
calculation (e), i.e. the one without mean field, gives completely different results 
for multiplicities and slope temperatures. This can also be seen in the light quark 
spectra given in Fig. |T^ where the results of calculations (a), (d) and (e) are shown 
by the solid, dashed and dotted hne, respectively. While (a) and (d) give almost 
identical results, calculation (e) differs substantially. Firstly, the energy range is 
shifted to smaller energies due to the small quark mass. Secondly, the spectrum 
is flatter, indicating a higher slope temperature. Thirdly, calculation (e) shows a 
deficiency of particles at low energies, instead of an enhancement like (a) and (d). 
This can be traced back to the fact that there is no mean field any longer, which 
shifts particles into this region. 

We have also investigated, whether the maximum in the temperature dependence 
of the relaxation times in Fig. |^ i. e. critical scattering, has some effect on the final 
spectra. To this purpose we have artificially set r(p, T) = r(p = 200 MeV, T) for 
p < 200 MeV and have compared the result with the full equation. The resulting 
spectra agree with those of Fig. |l^ within the accuracy of our calculation. 

We conclude that (i) medium effects in the collision terms have small effects on 
the shape of the final spectra but they considerably change the strangeness yield 
(here by 35%) and (ii) the medium effects in the mean field parts lead to significant 
effects in the spectra, especially in the low energy part. 
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IV. STRANGENESS PRODUCTION 



We discuss the variation of the strangeness yield with the radius tq of the initial 
fireball. The observables of the final state are the slope temperatures Tg and for 
the quarks and the strangeness yields. We have gradually increased the initial radius 
from To = 1 fm to ro = 7 fm while the other parameters of Eq. ( |3.31|) are kept fixed. 
The resulting values for the slope temperatures Tg and Tg and the strangeness yield 
Ms/Mq are given in Figs. ^ and |TB[ While the slope temperature for light quarks 
hardly changes with the system size, the slope temperature for the strange quarks 
rises and approaches Tq from below. 

The independence of Tq on tq corresponds to the independence of the hadronic 
slope temperatures on the projectile and target size, which has been reported in 
Ref. [0. The fact that T, stays below Tq seems to be in contradiction with ex- 
perimental data |jT9| and hydrodynamical models These models claim that 
the slope temperatures scale with the particle mass like T = Tq + m where 
(f^) is the mean transverse expansion velocity. One has to keep in mind, however, 
that this theoretical result is significantly influenced by the cylindrical geometry of 
the colliding nuclei and by the presence of transverse flow. These features are not 
present in our calculation. 

■ is essen- 



The increase of Tg with the initial radius tq, which can be seen in Fig. 15 



tially an effect of thermal equilibration. For small initial size, where there is negligi- 
ble thermal equilibration, one observes the "temperature" of the strange quarks di- 
rectly from the creation process. Its low value can be understood qualitatively by as- 
suming that the mean energy per particle E = mq + 3/2Tq = ms + 3/2Ts is conserved 
during the creation of a strange quark pair. One thus obtains Ts = Tq — 2/3{ms—mq) 
immediately after the pair creation. Subsequent elastic collisions, as they will occur 
more frequently in the larger systems, result in approaching Tq. This is indeed 
observed in Fig. 



According to Fig. 0, the strangeness yield grows with the radius of the initial 
flreball. This behaviour of the curve can be understood using the following simple 
analytical model. We make the ansatz 

^^'l^^' + div [vps{x, t)] = a [pg(f, t) - Kplix, t)] (4.1) 
for the space-time evolution of the strange quark density. In Eq. ( ^4.1| ) , the densities 
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are parametrized by 



P,iS^t) = J^Mu{t)-\x\) (4.2a) 



Ps{x,t) = ^^Mrs{t)-\x\) ■ (4.2b) 



For the radius of the hght quark system, we make the ansatz 

rl{t) = rl + vH' (4.3) 

with a constant velocity v. Equation ( [4.3|) is motivated by the exact solution of 
the Vlasov equation in the noninteracting case. We assume that the radius of the 
strange quark system behaves similarly, however, is smaller by a factor z: 

rs{t) = zr,{t) . (4.4) 

This factor accounts for the fact that only in the interior of the fireball the density 
is large enough to produce particles and is also born out numerically by comparing 
Figs. ^ and |^. With the parametrizations ( f4.2| )-( P^ ) at hand, Eq. ( |4.1| ) can be 
solved analytically to give 

. (4.) 

In the limit t ^ oo, one thus has 

§ . ,/Zta„h f!^) . (4.6) 
Since Mq ~ Tq, Eq. ( ^4.6| ) predicts a dependence on the initial radius as 

^ ~ tanh(aro) . (4.7) 

J\ q 



Figure shows, that this behaviour is indeed fulfilled in our calculation. The fit 
shown in Fig. |16| (solid line) has the parameters 

^ = 0.065 tanh(0.12ro) . (4.8) 

J\q 

The factor 6.5% in front of the hyperbolic tangent is the equilibrium value for an 
infinite system. A similar dependence like Eq. ( [4.7|) , however with different fit 
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coefficients, can be seen in the data for strangeness multiplicities at constant beam 



energy and variing mass number [^. The reason for the behaviour in Eq. ( [4 .71) is 
that for larger systems the interior zone keeps together longer and thus comes closer 
to chemical equilibrium. 

The asymptotic value of 6.5% in Eq. ( [4. 81 ) is about a factor of three small than 
the data given in Ref. There may be several reasons for this discrepancy. The 
ffist one is geometrical: Equation contains a factor z^/^, which accounts for the 
smaller size of the strange system compared to the light quark system. Estimating 
this factor to be 2; = 2/3 from our numerical results, this suppresses the strangeness 
yield by a factor two. It is probable that this factor will appear with a different power 
for a cylindrically symmetric system. The second reason for our low strangeness yield 
is that strangeness production in the initial high energy nucleon-nucleon collisions, 
like pp — > pAK, and strangeness production in the hadronic phase like irp KK 
are absent from our calculation. 



V. SUMMARY AND CONCLUSIONS 

We have developed numerical methods for the solution of transport equations 
whose physics is derived from the two and three flavor NJL Lagrangian and have 
given explicit solutions for various initial conditions. A finite difference algorithm is 
used to solve the partial differential equation for a spherically symmetric system. 

The paper is divided in two parts: First, the analytical and numerical properties 
of the two flavor Vlasov equation are given. Then this equation is generalized to 
include strangeness and also a collision term is added, so one deals with a Boltzmann 
type equation with selfconsistent medium dependent input. 

The most serious limitation encountered is the nonrenomalizability of the NJL 
model, which is treated by a cutoff of A = 602 MeV. For a thermal system with 
T = 240 MeV about 60% of all particles have momenta p > K and cannot be 
properly accounted for. 

For the Vlasov equation the total energy includes a potential as well as a kinetic 
term. The potential contributes ca. 5% to the total energy at t = 0. During the 
evolution, its contribution tends to zero as the mass evolves towards its vacuum 
value. The expansion dynamics in the Vlasov equation is dominated by the increase 
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of the constituent quark mass from rUq = 58 MeV at t = to = 312 MeV at 
t = oo. Correspondingly the quarks are slowed down during the expansion. For 
this reason the particle spectrum in momentum space changes during the expansion 
such that the low momentum region is enhanced and the high momentum region 
depleted. 

Simultaneously with the introduction of the collision term, the calculation is 
generalized to the three flavor model in order to study also the production of strange 
quarks. The collision term is implemented in a relaxation time approximation. 
In contrast to the Vlasov equation, the Boltzmann equation leads to an energy 
spectrum with an exponential shape at high energies. However, the mean field still 
manifests itself in an enhancement of the the low energy region of the spectra like 
for the Vlasov equation. The shapes of the spectra do not change when the medium 
effects in the cross sections are switched off. 

Only for large systems chemical as well as thermal equilibration of strange quarks 
can be reached. The approach to chemical equilibrium as a function of the initial 
radius tq can be well approximated by a tanh(aro) dependence. The absolute value 
of the strangeness yield found in our calculation is about a factor three smaller than 
observed in experiment [pT|. 

What have we learned from this calculation? Admittedly, an expanding quark 
system is an artificial system, since it lacks confinement. Yet it is a system in which 
the dynamics of a non-equilibrium system, which underlies chiral symmetry, can be 
studied. Medium effects are taken care of in the mean field as well as in the collision 
term. In this respect the sophistication surpasses other calculations using more 
realistic interaction models. While we cannot compare directly with experiment, yet 
at least three results may be of relevance: (i) The final spectra strongly depend on 
the mean field. For instance it is responsible for the enhancement at small energies. 
This is a direct consequence of the transition from a chirally symmetric to a chirally 
broken phase, (ii) The calculated spectra are rather robust with respect to changes 
in the collision terms. Whether one uses medium dependent cross sections or Born 
cross sections makes little difference, provided one does not modify the masses and 
the mean field, (iii) Medium effects are small for the masses of the strange quarks. 
The rather high values of rris are one reason, why strangeness production in the 
quark plasma alone accounts for only 30% of the observed value. 
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TABLES 

TABLE I. Particle multiplicities for different settings of initial system radius and max- 
imal momentum. In row (c), high momenta are cut off by a multiplicative factor. 









■!^q,t=oo 


J^s,t=oo 


i,t=oo/-^ j,t=oo 


{J^q+J^s)t=oo/Kt=0 


(a) 


1 


343.2 


365.6 


8.1 


2.2% 


1.089 


(b) 


3 


811.2 


825.4 


32.4 


3.9% 


1.058 


(c) 


1* 


251.2 


242.3 


5.7 


2.3% 


0.987 



TABLE IL Dependence of strangeness yield and slope temperatures on the medium 
effects. 





Mean Field 


Cross Sections 


A/'s,t=oo /A/'q,t=oo 


Tg (MeV) 


Ts (MeV) 


(a) 


yes 


exact 


2.2% 


185 


153 


(d) 


yes 


Born 


3.0% 


179 


153 


(e) 


no 


Born 


5.5% 


207 


191 
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FIGURES 

FIG. 1. The effective Potential Vefr as a function of the constituent quark mass niq. 
The parameters used for this plot are niQq = 5 MeV, A = 653.3 MeV and GA^ = 2.10. 
The constituent quark mass in the vacuum for this parameter set is 312 MeV. 

FIG. 2. The radial dependence of the quarks daensity as a function of r at various 
times t during the expansion process. 

FIG. 3. The radial dependence of the quark masses at various times t during the 
expansion process. 

FIG. 4. Momentum spectra for quarks in the beginning {t = 0, dashed line) and at 
the end (t = oo, solid line) of the expansion. 

FIG. 5. Time dependence of kinetic energy (dashed line) and total energy (solid line). 
Note the broken scale on the energy axis. 

FIG. 6. Inverse relaxation times (solid line), r^^^ (dashed line) and r^^^ (dotted 
line) as a function of temperature for p = 200 MeV. The dot-dashed line shows t~^^ in 
Born approximation. 

FIG. 7. Inverse relaxation time function of temperature for p = (solid line), 

p = 200 MeV (dashed line) and at p = in Born approximation (dotted line). 

FIG. 8. Light quark density at the same times as in Fig. ^ for the full calculation 
containing three flavors and collisions. 

FIG. 9. Strange quark density for the same times as in Fig. ^ for the full calculation 
containing three flavors and collisions. Note that the density is zero at t = and rises 
until t = A fm/c due to production by the process qq — > ss. 

FIG. 10. Light and strange quark masses as a function of time. 
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FIG. 11. Distributions of light and strange quarks for the final state as a function 
of the particle energy (solid lines). The dashed lines are exponential fits with the slope 
temperatures Tq = 185 MeV for light quarks and Tg = 164 MeV for strange quarks. 

FIG. 12. Distributions of light quarks resulting from different treatments of the high 
momentum part. Solid line: only particles with momentum p < A are included. Dashed 
line: particles with momentum p < 3A are included. Dotted line: smooth cutoff of high 
momentum particles. 

FIG. 13. Distributions of strange quarks resulting from different treatments of the high 
momentum part. Solid line: only particles with momentum p < A are included. Dashed 
line: particles with momentum p < 3A are included. Dotted line: smooth cutoff of high 
momentum particles. 

FIG. 14. Light quark spectra for three different treatments of medium effects; (a): 
standard procedure (solid line), (d): cross sections in Born approximation (dashed line) 
and (e): cross sections in Born approximation and mean field switched off. 

FIG. 15. Slope temperatures for light quarks (diamonds) and strange quarks (crosses) 
as a function of the radius ro initial fireball. 

FIG. 16. Variation of the strangeness yield as a function of the initial radius (dia- 
monds). The solid line is a fit by a hyperbolic tangent. 
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